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For a sound field observed on a sensor array, compressive sensing (CS) reconstructs the direction- 
of-arrival (DOA) of multiple sources using a sparsity constraint. The DOA estimation is posed 
as an underdetermined problem by expressing the acoustic pressure at each sensor as a phase- 
lagged superposition of source amplitudes at all hypothetical DOAs. Regularizing with an ^i-norm 
constraint renders the problem solvable with convex optimization, and promoting sparsity gives high- 
resolution DOA maps. Here, the sparse source distribution is derived using maximum a posteriori 
(MAP) estimates for both single and multiple snapshots. CS does not require inversion of the data 
covariance matrix and thus works well even for a single snapshot where it gives higher resolution than 
conventional beamforming. For multiple snapshots, CS outperforms conventional high-resolution 
methods, even with coherent arrivals and at low signal-to-noise ratio. The superior resolution of CS 
is demonstrated with vertical array data from the SWellEx96 experiment for coherent multi-paths. 


PACS numbers: 43.60.Pt, 43.60.Jn, 43.60.Fg 


I. INTRODUCTION 

Direction-of-arrival (DOA) estimation refers to the lo¬ 
calization of several sources from noisy measurements of 
the wavefield with an array of sensors. DOA estimation 
can be expressed as a linear underdetermined problem 
with a sparsity con stra int enforced on its solution. The 
compressive sensin^i^^ (CS) framework asserts that this 
is solved efhciently with a convex optimization procedure 
that promotes sparse solutions. 

In DOA estimation, CS achieves high-resolution acous¬ 
tic imagin^^l^^, outperforming traditional method^. 
Unlike the high-resolution subspace-based DOA 
estimatori^, DOA estimation via CS is reliable 
even with a single snapsholPtm. 

The least absolute shrinkage and selection operator 
(LASSO)P^ has been extended to mu ltiple measurement 
vectors (here multiple snapshots)PE^. They modify the 
LASSO objective function by introducing a mixed-norm 
penalty term that promotes spatial sparsity. More specif¬ 
ically, the snapshots are combined with the f’ 2 -norm, 
whereas the spatial samples are combined with the ii- 
norm. Multiple-snapshot CS offers s everal benefits over 
other high-resolution DOA estimator ^^111111 i) ft handles 
partially coherent arrivals. 2) It can be formulated with 
any number of snapshots, in contrast to, e.g., the Mini¬ 
mum Variance Distortion-free Response (MVDR) beam- 
former. 3) Its flexibility in formulation enables extensions 
to sequential processing, and online algorithmic^. Here, 
we show that CS achieves higher resolution than MUSIC 
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and MVDR, even in scenarios that favor these classical 
high-resolution methods. 

In ocean acoustics, CS has fo und several applications 
in matched field processin ^^^b^ l and in coherent passive 
fathometry for inferring sediment interfaces depths and 
their numbeilC^. Various wave propagation phenomena 
from a single source (refraction, diffraction, scattering, 
ducting, reflection) lead to multiple partially coherent ar¬ 
rivals received by the array. High-resolution beamformers 
cannot resolve these coherent arrivals. 


CS for single snapshot has high-resolution capabilities 
and contrary to eigenv alue-based beamformers works for 
coherent arrival^^^l. CS is limited by basis mismatchP^ 
which occurs when the DOAs do not coincide with the 
look directions of the angular spectrum, and by basis co¬ 
herence. Solutions to basis mismatch involve for exam¬ 
ple using atomic norm and solving the dual problenPE^ 
that are not addressed here. Grid refinement alleviates 
basis mismatch for high signal to noise ratio (SNR) at 
the expense of increased computational complexity. A 
denser grid causes increased coherence among the steer¬ 
ing vectors (basis coherence) which translates to bias 
and spread in the DOA estimates as demonstrated here. 
This is especially true in large two-dimensional or three- 
dimensional TOO-acoustic inversion problems as e.g. seis¬ 
mic imagingISHIl]. 


We use least squares optimization with an £i-norm reg¬ 
ularization term, also known as the LASSCS^, to formu¬ 
late the DOA estimation problem for single and multiple 
snapshots. The LASSO formulation complies with sta¬ 
tistical models as it provides a maximum a posteriori 
(MAP) estimate, assuming a Gaussian data likelihood 
and a Lapl acian prior distribution fo r the source acous¬ 
tic pressur^^mH fQj. both single (Sec. H.B) and multiple 
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snapshotJi^ (Sec. |lll[ ). The LASSO is known to be a con¬ 
vex minimization problem and solved efficiently by in¬ 
terior point methods. In the LASSO formulation, Sec. 
IV. A[ the reconstruction accuracy depends on the choice 


of the regularization parameter that controls the balance 
between the data fit and the sparsity of the solution. We 
indicate that the regularization para meter can be found 
from the properties of the LASSO pattP^I^, i.e., the evo¬ 
lution of the LASSO solution versus the regularization 
parameter. 

The main focus of the paper is on performance evalu¬ 
ation for single and multiple snapshots using both sim¬ 
ulated (Se c. |Vt and real data (Sec. VI). Other excel¬ 
lent papersld^ave already performed performance eval¬ 
uation for single snapshot, consistent with our simula¬ 
tions. We are not aware of performance evaluation for 
multiple snapshots. 

In the following, the £p-norm of a vector x G 

' ■ By extension. 


is defined as ||x|jp = 


the .^Q-norm is defined as ||x||o = '^^=1 and the 

£oo-norm as ||x||oo = max \xn\- For a matrix F G 

l<n<Ar 

f^MxL Frobenius norm H'lljr is defined as ||F||^ = 

E m I n |2 

i=l = ■ 


II. SINGLE SNAPSHOT DOA ESTIMATION 

We assume plane wave propagation and narrowband 
processing with a known sound speed. We consider 
the one-dimensional problem with a uniform linear array 
(ULA) of sensors with the source location characterized 
by the DOA of the associated plane wave, 0 G [—90°, 90°], 
with respect to the array axis. The propagation delay 
from the ith potential source to each of the M array sen¬ 
sors is described by the steering (or replica) vector. 


array signal-to-noise ratio (SNR) is defined as. 


SNR= lOlogio 


E{||Ax||i} 

E{||n||i} 


(dB). 


( 4 ) 


A. Sparse reconstruction with compressive sensing 

The problem of DOA estimation is to recover the set 
of non-zero components in the source vector x G C^, 
given the sensing matrix Amxn and an observation vec¬ 
tor y G . Even though there are only a few sources 
K < M generating the acoustic held, we are interested in 
a hne resolution on the angular grid to achieve precise lo¬ 
calization such that M N and the problem in Eq. (§ 
is underdetermined. A way to solve this ill-posed inverse 
problem is constraining the possible solutions with prior 
information. 

Traditional methods solve the problem in Eq. (U) by 
seeking the solution with the minimum £ 2 -norm which 
provides the best data ht (^ 2 -norm regularized least 
squares), 

Xf,(Ai) = argmin ||y - Ax\\l + fi\\x\\l. (5) 

xGC» 

The regularization parameter, /i > 0, controls the relative 
importance between the data ht and the ^ 2 -norm of the 
solution. The minimization problem in Eq. ([^ is convex 
with analytic solution, (AA^ -f /tIm) ^ y, 

where Im is the M x M identity matrix. However, it aims 
to minimize the energy of the source x through the £ 2 - 
norm regularization term rather than its sparsity, hence 
the resulting solution is non-sparse. 

Conventional beamforming (CBF)Pis related to the I 2 
solution for large /i. From Eq. (§: 

xcBF = lim inSie^in)) = A^y. (6) 

^—>■00 


a(6'i) = 


y/M 


1 , 


1 sin 9i 




, ( 1 ) 


where A is the wavelength and d the sensor spacing. 

Discretizing the half-space of interest, 6 G [—90°,90°], 
into N angular directions the DOA estimation problem 
can be expressed as a source reconstruction problem with 
the linear model. 


y = Ax -f n, (2) 

where y G is the complex-valued data vector from 
the measurements at the M sensors, x G is the un¬ 
known vector of the complex source amplitudes at all N 
directions on the angular grid of interest and n G is 
the additive noise vector. The sensing matrix, 

A= [a(0i),--- ,a(0^)], (3) 

maps the hypothetical sources x to the observations y 
and has as columns the steering vectors, Eq. Q, at all 
look directions. 

In the following, the noise is generated as independent 
and identically distributed (iid) complex Gaussian. The 


In principle, CBE combines the sensor outputs coherently 
to enhance the source signal at a specific look direction 
from the ubiquitous noise. CBF is robust to noise but 
suffers from low resolution and the presence of sidelobes. 

Since x is sparse (there are only K N sources), it 
is appropriate to seek for the solution with the minimum 
.^Q-aorm, which counts the number of non-zero entries in 
the vector, to find a sparse solution. However, the £ 0 - 
norm minimization problem is a non-convex combinato¬ 
rial problem which becomes computationally intractable 
ev en for moderate dimensions. The breakthrough of 
came with the proof that for sufficiently sparse 
signals, K N, and sensing matrices with sufficiently 
incoherent columns the 1 ' 0 -norm minimization problem 
is equivalent (at least in the noiseless case) to its con¬ 
vex relaxation, the £i-norm minimization probleirPfiUll 
By replacing the fp-norm with the convex £i-norm, the 
problem can be solved efficiently with convex optimiza¬ 
tion even for large dimensional^®!. 

For noisy measurements, Eq. ([^, the £i-norm mini¬ 
mization problem is formulated as 

Xf^(e) = argmin||x||i subject to jjy- Ax ||2 < e, (7) 

xGC™ 
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FIG. 1. Illustration of the LASSO path: Number of active 
indices versus the regularization parameter Increments in 
the active set occur at 


where e is the noise floor. The estimate (e) has the 
minimum ^i-norm while it fits the data up to the noise 
level. The problem in Eq. 0 can be equivalently written 
in an unconstrained form with the use of the regularizer 
/^ > 0 , 


Xfi(Ai) = argmin ||y - Ax||^ +/r||x||i. (8) 

xGC« 

The sparse source reconstruction problem in Eq. Q is a 
least squares optimization method regularized with the 
.^i-norm of the solution x and provides the best data fit 
(^ 2 -norm term) for the sparsity level determined by the 
regularization parameter /r. The optimization problem 
in Eq. Q is also known as the least absolute shrinkage 
and selection operator (LASSO) since the regularizer 
shrinks the coefficients of x towards zero as the regular¬ 
ization parameter fi increasei^^. This is illustrated in 
Fig.[g For every e there exists a /r so that the estimates 
in Eq. ([7) and Eq. ^ are equal. 

Once the active DOAs are recovered, by solving Eq. 0 
or equivalently Eq. (|^ , the unbiased complex source am¬ 
plitudes are determined from, 

xcs = A+^y, (9) 

where Am G contains only the “active” steering 

vectors associated with non-zero components in the solu¬ 
tion Xfj {fj.) and A'^ is its Moore-Penrose pseudoinverse. 

For a given sparsity level K and corresponding set of 
active indexes M, i.e. \M\ = A, Eq. © finds the best 
data fit. Thus, if the active sensing matrix A^vi has suf¬ 
ficiently incoherent columns it represents the solution to 
the £o problem 

X£q(A) = argmin ||y —Ax ||2 subject to ||x|jo = K (10) 

xGC« 


processes, by imposing a prior distrib ution on the solu¬ 
tion X which promotes sparsitj^^^^^^EH 

Bayes theorem connects the posterior distribution 
p(x|y) of the model parameters x conditioned on the 
data y, with the data likelihood p(y|x), the prior dis¬ 
tribution of the model parameters p(x) and the marginal 
distribution of the data p(y), 

Then, the maximum a posteriori (MAP) estimate is, 
xmap = argmax lnp(x|y) 

X 

= argmax [lnp(y|x) -|- Inp(x)] n2) 

= argmin [-lnp(y|x) - Inp(x)] , 


where the marginal distribution of the data p(y) is omit¬ 
ted since it is independent of the model x. 

Based on a complex Gaussian noise model with inde¬ 
pendent and identically distributed (iid) real and imag¬ 
inary parts, n ~ CA/'(0, cr^I), the likelihood of the 
data is also complex Gaussian distributed p(y|x) ^ 
CA/'(Ax, cr^I), 


P(y|x) = TT cr 


-Af„-2Afg- 


Ily-Axll 


(13) 


Followin^^, we assume that the coefficients of the solu¬ 
tion X are iid and follow a Laplacian-like distribution (for 
complex random variables). Such a prior has been shown 
to encourage sparsity in many situations because of the 
heavy tails and sharp peak at zero.The corresponding 
prior is 


N 

p(x) OC J|e 


y /(Re a:j)^+(Im xi)‘^ 


= e 


llxlli 


(14) 


The LASSO estimate, Eq. Q, can be interpreted as the 
MAP estimate. 


XMAP=argmin [||y - Ax\\l + fi\\x\\i]=Xi^{fi), (15) 

X 


where p. = jv. 

Equation (14) imposes no restriction on the source 
phases. Here, the phase is assumed uniformly [0, 27r) dis¬ 
tributed. 


III. MULTIPLE-SNAPSHOT DOA ESTIMATION 


B. MAP estimate via LASSO 

We use the LASSO formulation, Eq. to solve the 
DOA estimation problem in favor of sparse solutions. 
The choice of the (unconstrained) LASSO formulation 
over the constrained formulation, Eq. Q, allows the 
sparse reconstruction method to be interpreted in a sta¬ 
tistical Bayesian setting, where the unknowns x and the 
observations y are both treated as stochastic (random) 


Even though for moving sources it befits to solve one 
optimization problem for each snapshot sequentiallji^Sl^ 
for stationary scenarios, the sensor data statistics can be 
aggregated across snapshots to provide a more stable es¬ 
timate. Multiple snapshots are referred to as multiple 
measurement vectors and the recovery might have better 
performance than single measurement vectord^. Poten¬ 
tially the recovery can be made more robust by using 
a likelihood function Eq. (131 with colored noise (full 
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covariance matrix) or based on the Huber nornP^. For 
the multiple-snapshot case, all snapshot are collected into 
one matrix, 


Y = AX-bN, (16) 

where, for L snapshots, Y = [y(l),-- • ,y(i)] and N = 
[n(l), • • • , n(L)] are M x L matrices with the measure¬ 
ment and noise vectors per snapshot as columns, re¬ 
spectively, and X is the N x L signal with the com¬ 
plex source amplitudes at the N look directions per 
snapshot as columns. For stationary sources the matrix 
X = [x(l),- • • ,x(L)] exhibits row sparsity, i.e., it has a 
constant sparsity profile for every column, since the few 
existing sources are associated with the same DOA for all 
snapshots. As the sources are stationary it makes sense 
to sum the source energy across all snapshots, giving the 
row norm x^^ 

/L \l/2 


This quantity is sparse and in analogy with the single 
snapshot case we impose a Laplacian-like prior 


P(X) =p(x'^") (X exp(-||x^"||i/z/). 


(18) 


with no phase assumption. Similar to Eq. (14) we assume 
the phase is uniformly iid distributed on [0, 27r). 

We assume an iid complex Gaussian distribution for 
the data likelihood 


p(Y|X)cxexp(-||Y-AX||2./a2). (19) 

Using Bayes theorem, the MAP solution is then 
X = argmax p(Y|X)p(X) 

= argmin ||Y - AXjl^r-l-/r||x^^||i. (^*1) 


The beamformer power for CBF and MVDR respectively 
is then. 


Ucbf(6*) = wcbf(^)Cwcbf(^*) (23) 

^’mvdr(6*) = wmvdr(^')Cwmvdr(0), (24) 


where the corresponding weight vectors are given by. 


wcbf(^') = a(6>) 


wmvdr(0) = 


C-ia(6l) 


a^{e)c-^a{e)' 


(25) 

(26) 


The CBF can also be based directly on snapshots, as 
the single snapshot CBF Eq. ^ can be generalized to 
multiple snapshots, Xcbf = A^Y. The power estimates 
PcBpi^i), and the corresponding ith squared 

component of Xpg are thus comparable. Note that since 
the MVDR weights in Eq. (26) involve the inverse of the 


sample covariance matrix, MVDR requires a full rank C, 
i.e., L > M snapshots. 

The MUSICP is based on the eigendecomposition of 
the data sample covariance matrix Eq. (22) and the sep¬ 


aration of the signal and the noise subspace. 


C = U.A.U 


H 


U„A„U^. 


(27) 


The signal eigenvectors corresponding to the largest 
eigenvalues, Ag, are at the same subspace as the steer¬ 
ing vectors, Eq. Q, while the noise eigenvectors U„ are 
orthogonal to the subspace of the steering vectors thus 
a(0)'^U„ = 0. MUSIC uses the orthogonality between 
the signal and the noise subspace to locate the maxima 
in the spectrum. 


P 


MUSIC 


( 0 ) = 


a(0)^U„U^a(0)- 


(28) 


In this formulation we search for a sparse solution via 
the £i constraint. The source amplitude can, however, 
vary across snapshot. This is in contrast to covariance- 
matrix based beamforming that just inverts for the av¬ 
erage source power. The processing performance can be 
improved by doing an eigenvalue decomposition of X and 
retaining just the largest eigenvalues; see RefsEl IH The 
smaller eigenvalues contain mostly noise so this improves 
processing. However, this eigenvalue decomposition is 
not done here as this has features similar to forming a 
sample covariance matrix. 

Once the active steering vectors have been recovered, 
the unbiased source amplitudes are estimated for each 
snapshot, similar to the single snapshot case, Eq. (§, 

Xcs = AJY, (21) 

If desired, an average power estimate x^g can be obtained 
from the f 2 -iiorm of the rows of Xcs j with the ith element 
squared of x^g being the source power estimate at Oi. 

For reference, the CBF, MVDR, and MUSIC use the 
data sample covariance matrix, 

C = y YY^. (22) 

Lj 


Both MVDR and MUSIC overcome the resolution limit 
of the conventional beamformer by exploiting signal in¬ 
formation conveyed by the data sample matrix. How¬ 
ever, their performance depends on the eigenvalues of 
the data sample matrix thus it degrades with few snap¬ 
shots, when the data sample matrix is rank deficient, and 
in the presence of coherent sources, when the signal sub¬ 
space is reduced (Ch. 9 in Ref(7]). CS does not have these 
limitations as it utilizes directly the measured pressure 
Y. 


IV. REGULARIZATION PARAMETER SELECTION 

The choice of the regularization parameter /i in Eq. §, 
also called the LASSO shrinkage parameter, is crucial as 
it controls the balance between the sparsity of the esti¬ 
mated solution and the data ht determining the quality 
of the reconstruction. 

For large /r, the solution is very sparse (with small 
.^i-norm) but the data fit is poor as indicated in Fig. 

As /r decreases towards zero, the data fit is gradu¬ 
ally improved since the corresponding solutions become 
less sparse. Note that for ^ = 0 the solution in Eq. Q 
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becomes the unconstrained least squares solution. Since 
the LASSO path is derived and demonstrated for a single 
observation, the statistics of the source signal or noise is 
irrelevant. 


A. The LASSO path 


As the regularization parameter ^ evolves from oo to 
0, the LASSO solution in Eq. ^ changes continuously 
following a piecewise smooth tra jector y referred to as the 
solution path or the LASSO patlP^^^. In this section, we 
show that the singularity points in the LASSO path are 
associated with a change in the sparsity of the solution 
and can be used to indicate an adequate choice for /r. 

We obtain the full solution path using convex opti¬ 
mization to solve Eq. ([^ iteratively for different values 
of We use the CVX toolbox for disciplined convex op¬ 
timization that is available in the Matlab environment. 
It uses interior point solvers to obtain the global solution 
of a well-defined optimization probleirPSHSo] 

Let L(x,/r) denote the objective function in Eq. 

= l|y - Ax||^-h/r||x||i. (29) 


The value x minimizing Eq. (29) is found from its sub¬ 
derivative, 


d^L{x,fi) = 2A^ (Ax - y) -f ^5x||x||i, (30) 


where the subdifferential operator 9x is a generalization 
of the partial differential operator for functions that are 
not differentiable everywhere iReflSOl D.338j. The sub¬ 
gradient for the £i-norm is the set of vectors. 


9x||x||i = {s : ||s||oo < 1, s^x = |lx||i} , 
which implies. 


® “ |xi|’ 

|s^l < 1, 


Xiy^O 

Xi 0 , 


(31) 


(32) 


i.e., for every active element Xi ^ 0 oi the vector x G , 
the corresponding element of the subgradient is a unit 
vector in the direction of Xi. Eor every null element 
Xi = 0 the corresponding element of the subgradient has 
magnitude less than or equal to one. Thus, the magni¬ 
tude of the subgradient is uniformly bounded by unity, 

||s||oo < L 

Denote, 

r = 2A^ (y - Ax), (33) 


the beamformed residual vector for the estimated solu¬ 
tion X. Since Eq. (29) is convex, the global minimum is 
attained if 0 G /r) which leads to the necessary 

and sufficient condition 




G ^xllxlli. 


(34) 


Then, from Eq. (32) and Eq. (34), the coefficients = 
(y — Ax) of the beamformed residual vector r G 
have amplitude such that. 


\ri\=fi, Xi^O 
\ri\ < /r, Xi = 0, 


(35) 


i.e., whenever a component of x becomes non-zero, the 
corresponding element of the beamformed residual hits 
the boundary identified with the regularization parame¬ 
ter, ||r||oo < /r. 

For multiple snapshots, with the X determined from 
Eq. (20), the beamformed residuals become 


R = 2A^ (y - Ax) , 


Xi = 


\ 


El 

i=i 




(36) 


The values of fj, when changes in sparsity appear are ob¬ 
tained similarly to the single snapshot case. 


B. Algorithm for the LASSO path 


TABLE I. Fast iterative algorithm to solve the LASSO prob¬ 
lem (|^ for a desired sparsity level K and estimating the un¬ 
biased complex source amplitudes (§. 

Given: A G y G C^, A G N , E G]0, 1[ 

1: Initialize i = 0, = 0, = 2A^y 

2: while \Mi\ < K 
i = i -|- 1 

3: = (1 - F)peak(r‘-\ A) -b E peak(r*-\ A-b l) 

4a: x\^ — solution to Eq. ^ for A,y,y, = /r* 

4b = 2A^ {y-Ax\^) 

5. — {?n| l|co 

end 

6: if|Mi|>A 

7: Mi = {m| > 5i}, 5, = peak(|a;)J, A) 

end 

M=Mi 

8: xi^{y^) = x\ 

9: ®cs = A\^y 

10: Output: y'-, xcs, M. 


Although many algorithms exist for solving the LASSO 
problem, we have good experience with the algorithm in 
Table U as it is reasonable fast and accurate. Sec IIV.AI 
is used for formulating an algorithm where the values of 
y for different sparsity levels are indicated by the dual 
solution r, solving the dual problenP^. For large y, the 
solution X = 0 is trivial and r = 2A^y in Eq. (33). 
Decreasing y, a first component of x becomes active when 
the corresponding component of r hits the boundary, y = 
2|jA^y||oo, Xi = y. Inserting this solution into Eq. (33) 
and solving for the second peak of r hitting the boundary 
y indicates the value of y for which a second component 
becomes active. This way we follow the LASSO path 
in Fig. [l] towards less sparse solutions and lower y as 
detailed in Ref. [Ml 
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Starting from Eq. (331 with corresponding to 

regularization /i* for the set of active indexes M.i, the 
residual for the nth steering vector is now found. 

= 2a^ (y - Axe.ifi")) 


= 2a 


H 


y 


mGAii 


a. 




I y ^ ^ ) 

‘^a-nV 


(37) 

(38) 


The two progressively stronger approximations, in Eqs. 
([37|)-(|M|) , above are valid if the steering vectors corre¬ 


sponding to the fin al a ctive set is sufficiently incoherent 
la^Oml ~ 0. Eq. (38) actually corresponds to the con¬ 
ventional beamformer A^y for a single snapshot. The 
above equation is used for the selection of y, so it does 
not mean that the peaks in the conventional beamformer 
corresponds to the CS solution. 

The procedure is given in Table |lj where peak(r, fc) is 
the /cth peak of r. We choose F = 0.9. 

The dual method has been used to estimate the solu¬ 
tion path of the real-valuecP^ and the complex- valuecP^ 
generalized LASSO problems. The generalized LASSO 
uses the £i-norm to enforce structural or geometric con¬ 
straints on the solution by replacing the sparsity con¬ 
straint ||x||i with ||Dx||i for a structured matrix D. The 
generalized formulation performs well in certain applica¬ 
tions, e.g., recovery of continuous sources by promoting 
block sparsitji^ and DOA tracking for moving sources 
by an adaptive update of a diagonal weighting matrix 
D which reflects the evolution of the source probability 

distributiorPSl. 


C. Regularization parameter selection via the LASSO path 

The LASSO performance in DOA estimation is eval¬ 
uated by simulations starting with a large y and subse¬ 
quently decreasing its value. We consider an ULA with 
M = 20 sensors and spacing d = A/2. Three sources are 
at DOAs [—5, 0, 20]° with corresponding magnitudes [1, 
0.6, 0^ (linear) or [0,-4,—14] dB. The sensing matrix 
A in is defined on a coarse angular grid [—90°:5°:90°] 
(Fig. and a denser grid [—90°:1°:90°] (Fig. [^. The 
noise variance in (|^ is chosen such that SNR=20 dB. 

The trade-off between regularization term ||x||i and 
the data fit jjy — Ax||| in the LASSO estimate , Eq. (|^, 
for a range of values of y is depicted in Fig. The rele¬ 
vant values of y for the LASSO path are found between 
the two dots in Fig.j^b), i.e. 1.54 > y> 0.02. For these 
values of y, the importance shifts from favoring sparser 
solutions for large y towards diminishing the model resid¬ 
ual’s ^ 2 -norm for smaller y. From inspecting Fig. [^b), 
it is difficult inferring the value of y which results in the 
desired sparsity level. The LASSO path offers a more 
insightful method to determine the range of good val¬ 
ues of y (contained within the asterisks in Fig. [^b)) as 
explained below. 




FIG. 2. (Color online) The data error ||y — AxUl, describing 
the goodness of fit, versus the ^i-norm in (a) linear scale and 
(b) log-log scale for the estimated solution x for different val¬ 
ues of the regularization parameter y in the LASSO problem 
Eq. for sparse DOA estimation. 


Figure shows (a) the sparsity level ||x||o of the 
LASSO solution, (b) the properties of the LASSO path 
and (c) the corresponding residual vector versus the reg¬ 
ularization parameter y Since the interest is on sparse 
solutions X, it is natural inspecting the LASSO path for 
decreasing values of y, i.e., interpreting Fig. from right 
to left. 

For large values of y (e.g., y = 2) the problem in 
Eq. ^ is over-regularized, forcing the trivial solution 
X = 0 (Fig. I^b)), thus ||x||o = 0 (Fig. |3];a)). In this 
case, the slopes for all components jr-ij are zero (Fig.j^c)) 
since jr^j = |2af^y| < y for alH S [0, • • • , N] which is in¬ 
dependent of y. 

The first non-zero component of x appears at y = 
2||A^y|oo = 1-76 and remains active for y < 1.76 
(Fig. [3]b)) increasing the sparsity level to ||x||o = 
1 (Fig. [^a)). The corresponding component jr^j = 
\2aLf [y — anXi)\ (Fig. |^c)) is equal to y for ^ < 1.76. 
The other components Vj change slope at the singular 
point y = 1.76, since now |rj| = |2a^ (y — a^Xi)] < y 
for all j € [0, • • • ,A], j yf i. For y < 1.14, ||x|lo = 2 
(Fig. ^a)) as x acquires a second non-zero component 
(Fig. I^b)) and the corresponding component jr^j be¬ 
comes equal to y (Fig. [^c)). Similarly, the estimated 
solution has a third non-zero component for y < 0.38. 

For y < 0.18, x has many non-zero components 
(Figs.J^b),(c)) and its sparsity level increases abruptly 
(Fig. I^a)). For such low values of y the importance 
shifts to the data fitting term (^ 2 "iiorm term) in the reg¬ 
ularized problem, Eq. ([^, and x includes many non-zero 
noisy components gradually reducing the data error. 

The specific values of y at which an element of x be¬ 
comes active are denoted as the singular points in the 
piecewise smooth LASSO path. At a singular point, some 
component of r hits the boundary y, i.e. |r„| = y for 
some index n. Thus, the properties of the LASSO path 
indicate the selection of the regularization parameter y. 
For example, for a predefined sparsity level K a good 
choice of y is found by decreasing y until the Ath singu¬ 
lar point at the LASSO path. 

Owing to the piecewise smooth nature of the LASSO 
path, there is a range of y which give the same sparsity 
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FIG. 3. (Color online) The LASSO path versus /j, for three 
sources and SNR=20 dB. (a) Sparsity level of the estimate 
X. (b) Paths for each component of the solution x. (c) 
Paths for each component of the beamformed residual |r| = 
2|A^ (y — Ax) |. The vertical dashed lines indicates values 
of fi used in Figs. and 



FIG. 4. (Color online) The unbiased estimate xcs (°) for the 
true source x (*) and the corresponding beamformed residual 
vector for a coarse angular grid [—90:5:90]°. (a)-(b) /r = 

1.4, (c)-(d) /i = 0.5, (e)-(f) /i = 0.2, and (g)-(h) p = 0.1 
(corresponding to dashed lines in Fig.j^. The horizontal line 
in the residual plot (right) indicates the value of fi. 


level (i.e., between two singular points). In principle, the 
lowest fi in this range is desired as it gives the best data 
fit. Though, any value of fj, which achieves the desired 
sparsity suffices as once the active DOAs are recovered, 
the unbiased amplitudes are determined from Eq. 0- 

Figure 1^ shows the unbiased solution, Eq. along 
with the corresponding beamformed residual for four 
sparsity levels of /r. Notice how the residuals decrease 
in value as fj, is reduced. For ^ = 0.1, Fig. shows that 
five potential source locations exists as they have hit the 
boundary, so that components of |r| becomes equal to /r. 
Solving Eq. 0 shows that two sources are weak and are 
not shown. 

To increase precision in the LASSO reconstruction, a 
finer angular grid is required. However, angular grid 
refinement afso causes higher coherence among steering 
vectors, Eq. 0 : and the problem in Eq. 0 becomes 
increasingly underdetermined. Then, when solving the 
LASSO minimization Eq. ([^ might not exhibit the de¬ 
sired sparsity. Due to basis coherence and as /r decreases, 
components in the estimate x can be either activated 
(become non-zero) or annihilated. Similarly the resid¬ 
ual components can hit or leave the boundarj® (where 
components of |r| is equal to fi, see Eq. (HI)). 

In Fig. a denser angular grid with spacing 1° and 
setup as Fig. |^is used. For = 1.4, there is just one 
active component (Fig. [^a)) at —6° which is 1° away 
from the strongest DOA. This offset is mainly due to 
basis coherence, the correct location is not yet recovered. 
As /i is decreased the correct bin is eventually obtained 
(Fig-i e)). Thus when searching for a K sparse solution, 
it is often advantageous to search initially for more than 
K peaks, and then limit the final solution to the K most 
powerful elements. 
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FIG. 5. (Golor online) As in Fig.|^but for the denser angular 
grid [-90:1:90]°. 


The residual r is systematically decreased as /i is re¬ 
duced. All the active component can be seen as where 
components of jrj becomes equal to /i in the right panel 
of Fig. 


V. DOA ESTIMATION ERROR EVALUATION 

If the source DOAs are well separated with not too 
different magnitude, the DOA estimation for multiple 
sources using CBF and CS turns out to behave simi¬ 
larly. They differ, however, in their behavior whenever 
two sources are closely spaced. The same applies for 
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MVDR under the additional assumptions of incoherent 
arrivals and sufficient number of snapshots, L > M. The 
details are of course scenario dependent. 

For the purpose of a quantitative performance eval¬ 
uation with synthetic data, the estimated, 9^, and the 
true, DOAs are paired with each other such that 

the root mean squared DOA error is minimized in each 
single realization. After this pairing, the ensemble root- 
mean-squared error is computed. 


RMSE = 




E 


1 ^ 
K ^ 

fc=i 


(0fc - 


(39) 


The data is generated to have a fixed SNR Eq. 0. 
The source phases of each x component is uniformly dis¬ 
tributed on [0, 2tt) in order to generate a sample covari¬ 
ance matrix from which MVDR/MUSIC can resolve in¬ 
coherent sources. 

CBF suffers from low-resolution and the effect of side- 
lobes for both single and multiple data snapshots, thus 
the simple peak search used here is too simple. These 
problems are reduced in MVDR and MUSIC for multiple 
snapshots and they do not arise with CS. 

The optimal performance for K sources is found by 
searching over all combinations of steering vectors for the 
maximum likelihood solution, Eq. (131, i.e., the best fit¬ 
ting source vector using Eq. (|^. This is a NP-hard com¬ 
binatorial problem, that for N look directions requires 
evaluation of N\/K\{N — K)\ solutions. For N=361 and 
either K=2 or K=3, this gives 77,000 or 7,700,000 com¬ 
binations to be evaluated. This makes the exhaustive 
search approach impractical for larger K. 

In the following simulation, we consider an array with 
M = 20 elements with spacing d = A/2. The DOAs are 
assumed to be on a fine angular grid [—90°:0.5°:90°], i.e. 
A G ^20x361^ regularization parameter /i is chosen 
to correspond to the K + 2 largest peak of the residual 
in Eq. (33) using the procedure in Table |T] and retaining 
only the K largest source powers. We require the peaks 
of the CS to be at least 4 bins apart. Thus the exhaustive 
and the CS do not solve the identical problem, as the CS 
solves a smaller problem. Note that panel c in Figs. |6}|^ 
shows the simulation results versus array SNR defined in 
Eq. (||). 


A. Single Snapshot 
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FIG. 6. (Color online) Single snapshot example for 2 sources 
at DOAs [2, 75]° and magnitudes [22, 20] dB. At SNR = 5dB 
a) spectra for CBF, CS (o) and unbiased CS (o, higher 
levels), and b) CS, CBF and exhaustive-search histogram 
based on 1000 Monte Carlo simulations, and c) CS, CBF and 
exhaustive-search performance versus SNR. The true source 
positions (*) are indicated in a) and b). 


are characterized by a high sidelobe level but for the two 
well-separated similar-magnitude sources this is a minor 
problem. 

Using Monte Carlo simulations, we repeat the CS in¬ 
version for 1000 realizations of the noise in Fig. The 
RMSE increases towards the endfire directions. This is 
to be expected as the main beam becomes wider and this 
results in a lower DOA resolutiorP. Since the sources are 
well-separated in this scenario, CS, CBF, and exhaustive 
search perform similarly with respect to RMSE. 


In the first scenario, we consider a single snapshot case 
with additive noise with K = 2 well-separated DOAs at 
[2, 75]° with magnitudes [22, 20] dB, see Fig. In the 
second scenario, a third weak source is included very close 
to the first source: Thus, K = 3 and the source DOAs 
are [—3, 2, 75]° with magnitudes [12, 22, 20] dB, see Fig. 

The synthetic data is generated according to Eq. ([^. 

For the first scenario, the CS diagrams in Fig. 0 show 
DOA estimation with small variance but indicate a bias 
towards endfire, as for the true DOA 75° the CS estimate 
is 76°. Towards endfire the main beam becomes broader 
and absorbs more noise power. The CBF spectra Fig. §1 


Repeating the Monte Carlo simulations at several 
SNRs gives the RMSE performance of CS and CBF in 
Fig. [^. Their performance is about the same since the 
DOAs are well-separated. 

In the second scenario, the CBF cannot resolve the 
two closely spaced sources with DOAs [—3, 2]°. They 
are less than a beamwidth apart as indicated in Fig. 0- 
Sidelobes cause a few DOA estimation errors at —65° in 
the CBF histogram. Fig. [^. Since CS obtains high- 
resolution even for a single snapshot, it performs much 
better than CBF, Fig. [^. 
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FIG. 7. (Color online) As Fig. but for 3 sources at DOAs 
[—3, 2, 75]° and magnitudes [12, 22, 20] dB. 


B. Multiple Snapshot 


In the multiple-snapshot scenario, MVDR and MU¬ 


SIC use the data sample covariance matrix Eq. (22) 


whereas CBF and CS works directly on the observations 
X Eq. (201. The sample covariance matrix is formed by 


averaging L synthetic data snapshots. The source magni¬ 
tude is considered invariant across snapshots. The source 
phase is for each snapshot sampled from a uniform dis¬ 
tribution on [0, 27r). 

Due to the weak performance of MVDR in scenarios 
with coherent arrival^, we assume incoherent arrivals 
in the simulations although not needed for CS. For CS 
we use Eq. (20) with a similar choice of regularization 


parameter as for the single snapshot case. 

Using the same setup as in Fig. but estimating the 
source DOAs based on L = 50 snapshots gives the results 
in Fig. At SNR = 0 dB the diagrams in Fig. show 
that CS localizes the sources well, in contrast to the CBF 
and MVDR that is also indicated in the histograms in 
Fig. [^. The RMSE in Fig. [^, shows that CBF does not 
give the required resolution even for high SNR. MVDR 
performs well for SNR > 10 dB, whereas CS performs 
well for SNRs down to 2.5 dB. 

In a third scenario, the weak broadside sources are 
moved closer with DOAs defined as [—2, 1, 75]°. Fig. 
gives about the same DOA estimates for CBF, as it is 
already at its maximum performance even for high ar¬ 
ray SNR, confirming its low resolution. MVDR fails for 
SNR < 20 dB, which is 10 dB higher than the correspond¬ 


ing value in Fig. (MUSIC fails also at a level 10 dB 
higher). Contrarily, CS fails only for SNR < 5 dB which 
is 2.5 dB higher (Figs.[^ and|^). Note how MVDR com- 
pletely misses the weak source at —2° in Figs. [^, but 
CS localize it with a larger spread. Thus, as the weak 
source moves closer to the strong source, CS degrades 
slower than MVDR in terms of RMSE. This is a good 
indication of its high-resolution capabilities. 

Figure IT shows the estimated power at the one re¬ 
alization in Fig. [§1 of L = 50 snapshots inverted si¬ 
multaneously. We emphasize the scale of the problem. 
Equation ( [T^ has 20 • 50 = 1000 equations to determine 
361-50 = 18050 complex-valued variables at 361 azimuths 
and 50 snapshots observed on 20 sensors. The sparsity 
constraint is crucial here. 


The CS (and especially the exhaustive-search) requires 
several orders of magnitude more CPU-time than the 
beamforming methods. 

Many other simulations could be performed, for exam¬ 
ple colored noise, no assumptions on number of sources, 
and random source locations. From initial exploration 
of these it is our impression that CS will perform well, 
though more simulations are required. 
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FIG. 8. (Color online) Multiple L = 50 snapshot example for 
3 sources at DOAs [—3, 2, 75]° with magnitudes [12, 22, 20] 
dB. At SNR = OdB a) spectra for CBF, MVDR, and CS 
(o) and unbiased CS (o, higher levels), and b) CS, CBF and 
MVDR histogram based on 100 Monte Carlo simulations, and 
c) CS, exhaustive-search, CBF, MVDR, and MUSIC perfor¬ 
mance versus SNR. The true source positions (*) are indicated 
in a) and b). 
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FIG. 9. (Color online) As Fig. but with closer spaced 
sources [—2, 1, 75]°. 



FIG. 10. (Color online) Power (linear) for the multiple 
snapshot case across azimuths and snapshots for one noise 
realization at SNR = OdB for the scenario with DOAs at 
[-3, 2, 75]°. 


VI. EXPERIMENTAL RESULTS 

The high-resolution performance of CS both in single- 
and multiple-snapshot cases is validated with experimen¬ 
tal data in a complex multi-path shallow-water environ¬ 
ment and it is compared with conventional methods, 
namely CBF and MVDR. 

The data set is from the shallow wat er eva luation cell 
experiment 1996 (SWellEx-96) Event collected on 

a 64-element vertical linear array. The array has uni¬ 
form intersensor spacing 1.875 m and was deployed at 
waterdepth 16.5 m spanning 94.125-212.25 m. During 
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FIG. 11. (Color online) Spatial CBF spectrum across fre¬ 
quency at the source’s closest point of approach to the array. 


the Event S5, from 23:15-00:30 on 10-11 May 1996 west 
of Point Loma, CA, two sources, a shallow and a deep, 
were towed simultaneously from 9 km southwest to 3 km 
northeast of the array at a speed of 5 knots (2.5 m/s). 
Each source was transmitting a unique set of tones. 

Here, we are interested in the deep source towed at 
54 m depth while at the vicinity of the closest point of 
approach (CPA) which was 900 m from the array and 
occurred around 00:15, 60 min into the event. The deep- 
towed source signal submitted a set of 9 frequencies [112, 
130, 148, 166, 201, 235, 283, 338, 388] Hz at approxi¬ 
mately 158 dB re 1/rPa. The processed recording has 
duration of 1.5 min (covering 0.5 min before and 1 min 
after the CPA) sampled at 1500 Hz. It was split into 87 
snapshots of 2^^ samples (2.7 s) duration, i.e., with 63% 
overlap. 


Figure [m shows the multiple-snapshot CBF spatial 
spectrum, Eq. (24), over the 50-400 Hz frequency range. 
Arrivals are detected not only at the transmitted tonal 
frequencies of the deep towed source but also at sev¬ 
eral other frequencies corresponding to the shallow-towed 
source tonal frequencies, weaker deep source frequencies, 
and the acoustic signature of the tow-ship. 

Single-snapshot processing with CBF and CS at the 
deep source tonal set, contour plots in Fig. [T^ indicates 
the presence of several multipath arrivals which are ade¬ 
quately stationary along the snapshots at the CPA. Due 
to the significant sound speed variation it is not straight¬ 
forward to associate the reconstructed DOAs with spe¬ 
cific reflections. The CBF map comprises 6 signihcant 
peaks but suffers from low resolution and artifacts due 
to sidelobes and noise. To choose the regularization pa¬ 
rameter in the LASSO formulation for CS reconstruction, 
we solve iteratively Eq. (j^ as described in Table jljwith 
initial value /r = 2])A^y|[^, until the obtained estimate 
has a sparsity level of 10. The CS reconstruction results 
in improved resolution due to the sparsity constraint and 
significant reduction of artifacts in the map. 

Combining the data from all the snapshots and pro¬ 
cessing with CBF, MVDR, and CS, as in Sec |V.B[ re¬ 
veals that MVDR fails to detect the coherent multipath 
arrivals; see line plots in Fig.[^ Again the peaks of CBF 
and CS are consistent but CS offers improved resolution. 

We have here used higher sparsity for the single¬ 
snapshot processing to allow for identifying non- 
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FIG. 12. (Color online) Single (contour plots) and multiple (line plots) snapshot reconstruction at the transmitted frequencies 
with CS (★), CBF (background color, solid) and MVDR (dashed). For the single snapshot we have assumed 7F = 10 sources 
while for the multiple snapshot K = 6. 


stationary paths. The non-stationary path can be seen 
in several of the contour plots, most prominently at 112, 
130 and 201Hz. When performing multiple-snapshot pro¬ 
cessing where the solution is constrained to remain active 
at one azimuth (but with varying power), the stationary 
paths are most likely to contribute to the CS solution. 


VII. CONCLUSION 

The estimation of multiple directions-of-arrival (DOA) 
is formulated as a sparse source reconstruction problem. 
This is efficiently solvable using compressive sensing (CS) 
as a least squares problem regularized with a sparsity 
promoting constraint. The resulting solution is the max¬ 
imum a posteriori (MAP) estimate for both the single 
and multiple-snapshot formulations. The regularization 
parameter balances the data fit and the solution’s spar¬ 
sity. It is selected so that the solution is sufficiently sparse 
providing high-resolution DOA estimates. A procedure 
to find an adequate choice for the regularization param¬ 
eter is described whereby the DOAs are obtained. 

CS provides high-resolution acoustic imaging both 
with single and multiple snapshot. The performance 
evaluation shows that for single snapshot data, CS gives 
higher resolution than CBF. For multiple snapshots, CS 
provides higher resolution than MVDR/MUSIC and the 
relative performance improves as the source DOAs move 
closer together. 

The real data example indicates that CS is capable of 
resolving multiple coherent wave arrivals, e.g. stemming 
from multipath propagation. 
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